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Abstract —Approximate message passing (AMP) and its vari¬ 
ants, developed based on loopy belief propagation, are attractive 
for estimating a vector x from a noisy version of z = Ax, which 
arises in many applications. For a large A with i. i. d. elements, 
AMP can be characterized by the state evolution and exhibits 
fast convergence. However, it has been shown that, AMP may 
easily diverge for a generic A. In this work, we develop a new 
variant of AMP based on a unitary transformation of the original 
model (hence the variant is called UT-AMP), where the unitary 
matrix is available for any matrix A, e.g., the conjugate transpose 
of the left singular matrix of A, or a normalized DFT (discrete 
Fourier transform) matrix for any circulant A. We prove that, in 
the case of Gaussian priors, UT-AMP always converges for any 
matrix A. It is observed that UT-AMP is much more robust than 
the original AMP for ‘difficult’ A and exhibits fast convergence. 

A special form of UT-AMP with a circulant A was used in our 
previous work fl3l for turbo equalization. This work extends it 
to a generic A, and provides a theoretical investigation on the 
convergence. 

Index Terms —Belief propagation, approximate message pass¬ 
ing (AMP), convergence, singular value decomposition (SVD). 

I. Introduction 

PROXIMATE message passing, developed based on 
loopy belief propagation, is an efficient approach to the 
estimation of a vector x with independent elements {xi ~ 
p(xi)} in the following model 

y = Ax + n (1) 

where A is a known matrix with size M x N, the length of 
x is N, y denotes a length-M observation vector and n is a 
length-M white Gaussian noise vector with zero mean and 
covariance matrix er 2 I 01-0. AMP was originally developed 
for compressive sensing based on model 01 ffl-l3l, and then 
was extended to generalized AMP (GAMP) to accommodate 
more general distribution p(j/j|(Ax)j) which may not be Gaus¬ 
sian (where and (Ax), denotes the z-th element in y and 
(Ax), respectively) (4), 0. For a large A with i.i.d. elements, 
AMP exhibits fast convergence which can be characterized 
by the state evolution 0, 0. However, for a generic A, the 
convergence of AMP cannot be guaranteed. It has been shown 
that AMP may diverge for a benign matrix A and it can easily 
diverge for a ‘difficult’ matrix A, e.g., non-zero mean, rank- 
deficient, column-correlated, or ill-conditioned A 0, ED- 

The fixed points and convergence of AMP were analyzed for 
an arbitrary matrix A in 0 and EQ . Reference 10 provides 
sufficient conditions for the convergence of AMP in the case of 
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Gaussian priors {p(xi)}. The convergence condition is closely 
related to the peak-to-average ratio of the squared singular 
values of a certain normalized A for vector stepsize AMP 
algorithm, and is closely related to the peak-to-average ratio 
of the squared singular values of A for scalar stepsize AMP 
algorithm. Damped AMP algorithms were proposed and the 
convergence can be guaranteed with sufficient damping, but 
the amount of damping grows with the peak-to-average ratio 
0. Adaptive damping and mean removal mechanisms were 
introduced to (G)AMP in (T2) to enhance the convergence 
speed. Compared to original AMP, swept AMP (SwAMP) in 
0, ED is much more robust to difficult A. However, SwAMP 
updates the relevant estimates sequentially (in contrast, AMP 
updates them in parallel), which restricts fast implementations. 
The global convergence of AMP with a generic A for a generic 
prior p(xi) has not been understood 0. 

In this work, we present a new variant of AMP, which is 
developed based on the following unitary transformation of 
0 

r = AV + w (2) 

where r = U ff y, w = U^n, and 

A = UAV (3) 

with U and V being unitary matrices and A being a rectangular 
diagonal matrix. We note that, as T5 H is a unitary matrix, 
w is still a zero mean Gaussian noise vector with the same 
covariance matrix as n in 0. Eqn. 0 holds for any A 
through the singular value decomposition (SVD). It is worth 
mentioning that, any circulant matrix A (M = N) can be 
unitarily diagonalized by a discrete Fourier transform (DFT) 
matrix, so U and V can simply be the normalized DFT matrix 
and its inverse. In addition, r and (the diagonal elements of) 
A can be calculated with the fast Fourier transform (FFT). A 
new variant of AMP is then developed based on model 0, 
which, for convenience, is called UT-AMP (where UT stands 
for unitary transformation) in this paper. It is interesting that, 
although unitary transformation does not change the singular 
values of A, we prove that UT-AMP converges for any A 
in the case of Gaussian priors. Moreover, we show that the 
convergence speed of UT-AMP is related to a scalar a (see 
Theorem 1 and its proof). It is observed that UT-AMP is much 
more robust than the original AMP algorithms and exhibits 
fast convergence. It is noted that the SVD required for a 
non-circulant A only needs to be carried out once, so UT- 
AMP is particularly suitable for applications with a fixed A 
(e.g., turbo MIMO detection with/without quasi-static channels 
in communications). For applications with model 0 with 
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a circulant A (e.g., block transmission with cyclic prefix in 
communications), the unitary transformation can be efficiently 
performed with FFT, which makes UT-AMP very attractive, 
e.g., in equalization to combat intersymbol interference, as 
shown in our previous work fj~3l . 

Notations: Bold lowercase letters, e.g., c, are used to denote 
column vectors, and bold upper case letters, e.g., C, are used 
to denote matrices. The i-th element in vector c is denoted by 
d. We use c • d and c./d to denote the elementwise product 
and division between two vectors c and d, respectively. |C| 2 3 4 5 6 * 8 9 
represents the elementwise magnitude squared operation for 
matrix C. 1, 0 and I represent an all-one column vector, an all¬ 
zero column vector, and an identity matrix with proper sizes 
depending on the context. The conjugate transpose is denoted 
by the superscript “ H ”. 


II. AMP WITH VECTOR STEPSIZES AND SCALAR 
STEPSIZES 


To facilitate comparisons with UT-AMP, we include the 
vector stepsize AMP (Algorithm 1) and the scalar stepsize 
AMP (Algorithm 2) 0 in this section. In vector stepsize 
AMP, the function g x (q. T q ) returns a column vector whose 
i -th element, denoted by [g x (q.T q )]i, is given by 


[dx (<T T q)\i 


f Xip{xi)N[xj] quTg^dxj 
f p(xi )Af{xi ; qi , r qi )dxi 


(4) 


where Af(xi] qi,T qi ) denotes a Gaussian distribution with Xi 
as random variable, q, as mean, and t, h as variance. Eqn. © 
can be interpreted as the MMSE (minimum mean square error) 
estimation of Xi based on the following model 


qi = Xi+u7 (5) 

where a \ ~ p(xi) and vj is a Gaussian noise with mean 
zero and variance r qi . The function g' r (q . T q ) returns a column 
vector, and the i-th element is denoted by [g x (q, T q )\i where 
the derivative is with respect to qi. It is not hard to show that 
T qi \g' x { q, T q )\i is the a posteriori variance of Xi with model ©. 
Note that g x ( q, T q ) can also be changed for MAP (maximum 
a posteriori) estimation of x 0. 

Scalar stepsize AMP can be obtained from vector stepsize 
AMP by forcing the elements of each variance vector to be the 
same, so that the multiplications of a matrix with a vector in 
updating t p and r q are avoided (compare Lines 1 and 5 in both 
algorithms). The function g x ( q, r q ) is the same as g x (q, T q ) in 
vector stepsize AMP except that all the Gaussian distributions 
{Af(xi‘, qi, r q )} share the same variance r q . 


III. UT-AMP and Its Convergence 
A. Derivation of UT-AMP 

As any matrix A can have the decomposition A = UAV, we 
first perform a unitary transformation with \J H to ©, yielding 

r = U H y = (AV)x + w (6) 

where A is an M x N rectangular diagonal matrix. Then the 
vector stepsize AMP can be applied to © where the system 
matrix becomes a special matrix AV. Note that 

|C| 2 d = (C Diag{ d) C H ) D 1 (7) 


Algorithm 1 Vector Stepsize AMP 

Initialize tJ°' ) (with elements larger than 0) and x-'U Set 
s( -1 ) = 0 and t = 0. 

Repeat 

1. T P = |A| 2 t* 

2. p = Ax* — T p • s* _1 

3. T s = 1 ./{Tp + cr 2 l) 

4. s* = T, ■ (y - p) 

5. l./Tq = |A ff | 2 r s 

6. q = x* + T q ■ A H s t 

7 ■ T x +1 = T q -g' x (q,T q ) 

8. x t+1 = g x (q,T q ) 

9. t = t+l 

Until terminated 


Algorithm 2 Scalar Stepsize AMP 

Initialize ri °' * 1 > 0 and x^°\ Set = 0 and t = 0. 

Repeat 


1. Tp = 

(1/M)|A| 2 f t* 

2 P = 

Ax* — TpS* 1 

II 

l/(Tp + (J 2 ) 

4. s* = 

^(y-p) 

5- 1 /T ? 

= (1/N)\A h \ 2 f t s 

6 . q = 

X* + Tq A ff S* 

7 - 4 +1 

= {r q /N)l H g' x (q,T g ) 

8 . x * +1 

= 9x{q,T q ) 

9. t = t+ 1 

Until terminated 


where Diag(d) returns a diagonal matrix with the elements 
of d on its diagonal, and ( B ) /; returns a diagonal matrix by 
forcing the off-diagonal elements of B to zero. Now suppose 
we have a variance vector r*. According to Line 1 in vector 
stepsize AMP and using ©, we have 

T p = (AV Diag{T x ) \ H A h )t>1. (8) 

In attempting to reduce the computational complexity, we can 
find that if rj, has a form of 7 I, the calculation of © can be 
significantly reduced. This motivates the replacement of rj. 


Algorithm 3 UT-AMP 

Unitary transform: r = U^v = AVx + w, where A = UAV. 
Define vectors A p = AA H 1 and A s = A H Al. 

Initialize r|°^ > 0 and x^°\ Set s^ -1 ^ = 0 and t = 0. 

Repeat 

1. Tp = T x Xp 

2. p = AVx* - t p ■ s*’ 1 

3 . T s = 1 ./(Tp + cr 2 l) 

4. s* = T s ■ (r - p) 

5. 1 /Tq = (1/JV)A? T S 

6. q = x* + r q (Y H A H s t ) 

7 ■ ri +1 = {T q /N)l H g' x (q,r q ) 

8. x t+1 = g x (q. T q ) 

9. t = t+ 1 
Until terminated 
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with r* 1 where r* is the average of the elements of t,*. So 
([S j is reduced to 

r.p = r‘AA ff l (9) 


B. Convergence of UT-AMP 

Theorem 1. UT-AMP converges for any A in the case of 
Gaussian priors. 


which is Line 1 in UT-AMP. Lines 2, 3 and 4 in UT-AMP can 
be obtained according to Lines 2, 3, 4 in vector stepsize AMP 
by simply replacing A with AV. According to (|7]i again. Line 
5 in vector stepsize AMP with matrix AV can be represented 
as 

l./rq = (Y H A H Diag(r p ) AY) D 1. (10) 

In order to reduce the computational complexity, we can 
replace the diagonal matrix A H Diag(r p ) A in ( fTTTb with a 
scaled identity matrix pi where p is the average of the diagonal 
elements of A H Diag(r p ) A, i.e., 

P =(1/N)1 h A h At p . (11) 

Hence (ITOt is reduced to Line 5 in UT-AMP. Line 6 can be 
obtained from Line 6 in vector stepsize AMP by replacing A 
with AV. Compared with Line 7 in vector stepsize AMP, an 
additional average operation is performed in Line 7 in UT- 
AMP to meet the requirement of a scalar r* in Line 1. We 
note that the average operation is not necessarily in Line 7 as 
we can also put the additional average operation in Line 1. 
Line 8 in UT-AMP is the same as Line 8 in vector stepsize 
AMP except that r q is a scalar. 

Remarks: 

• One may try to get another variant of AMP by applying 
the scalar stepsize AMP to model i.e., replacing A 
with \} H A and replacing y with r = U H y in scalar step- 
size AMP. It is interesting that the obtained algorithm will 
remain exactly the same as the original scalar stepsize 
AMP as U /; will be canceled out in scalar stepsize AMP. 
This means that unitary transformation have no impact on 
the convergence of scalar stepsize AMP. 

• By the name, in vector stepsize AMP, r*, r p , r s , and 
r q are all vectors, and in scalar step size AMP, the 
corresponding r*, t p , t s , and r q are all scalars. In 
contrast, UT-AMP has two scalars r* and r q and two 
vectors r p and t s . 

• If A is a circulant matrix, UT-AMP is very attractive as 
U and V can be simply a DFT matrix and its inverse, and 
the diagonal elements of A can be calculated with FFT. 
Moreover, the multiplications of matrix and vector in UT- 
AMP can be implemented with FFT as well, leading to 
very low complexity. 

• If A is non-circulant and its SVD required in UT-AMP 
is available, the complexity per iteration of the UT- 
AMP is lower than that of vector stepsize AMP as 
the multiplications of matrix with vector are avoided in 
Lines 1 and 5. The complexity of UT-AMP is slightly 
higher than that of the scalar stepsize AMP due to the 
vector operations in Lines 1 and 5. Hence, UT-AMP is 
particularly suitable for applications with fixed A as SVD 
only needs to be carried out once. 

• Most importantly, it is observed that UT-AMP is robust 
to ‘difficult’ matrix A and exhibits fast convergence. 


Proof: See Appendix A. ■ 

It can be seen from the proof of Theorem 1 that, the 
convergence speed of UT-AMP is related to a parameter a 
given in (l30l) in Appendix A. 

Similar to the original AMP, the convergence of UT-AMP 
for a generic prior is unknown, which remains as future work. 
It is also interesting to investigate the convergence of swept 
UT-AMP. 

It is observed that UT-AMP is robust to ‘difficult’ A 
e.g., non-zero mean, rank-deficient, column-correlated, or ill- 
conditioned A, under which the original AMP often diverges. 
The special form of UT-AMP with a circulant A in the 
case of discrete priori distributions has been used in lfl3l for 
equalization, where the channel matrix A is ill conditioned, 
and the use of original AMP will diverge. Various numerical 
examples will be provided in the full version of this paper. 

IV. Conclusion 

In this work, we have developed a new AMP variant UT- 
AMP for a generic matrix A. It has been shown that UT-AMP 
always converges in the case of Gaussian priors for any A. It 
is observed that UT-AMP is robust to difficult A and exhibits 
fast convergence. 


Appendix A 
Proof of Theorem 1 

We assume that 

p(x) ~ Af(\-,\ 0 ,Diag{r°)) (12) 

where x° and r° are the a priori mean vector and variance 
vector for x, respectively. 

Similar to the proof in || 8 ], it can be proven that the variance 
t* of UT-AMP for any A always converges to a fixed point 
denoted by t x . Next, we prove the convergence of x t . 

Define a diagonal matrix 

V = Diag{T s ) = {tIAA h + a 2 !)- 1 . (13) 

Then, according to the UT-AMP algorithm, we have 

s* = t*DAA h s‘ -1 - DAVx* + Dr, (14) 


l/ Tq = {1/N)1 h A h BA1 


min{M,N } 

= 1 y 

N ^ 


I Ail 


r| M 2 + a 2 


with A i being the (*,i)-th elements of A. and 


yt+l 


x” 1 " = r* +1 (q/r 9 +x°./r“) 

= (r‘ + 1 /t 9 )x‘ + r t + 1 Y H A H S t + t* + 1 x°./t° 


(15) 


= T, 


t+1 rlY H A ff DAA ff s t_1 


(cil - r* +1 V ff A ff DAV)x t + b, (16) 
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where b is an appropriate vector. Define 

_ r; 


a = 


.4+1 I t+ i|. i 2 

a: _ 1 ^ T x l A »l 


r 9 iV ^ T*|Aj| 2 + cr 2 


(17) 


The iteration of x* and s* in UT-AMP can be described as 


s* 


' Ca 

c 6 ' 


s* 1 

x t+1 


. C c 

c d 


X* 


(18) 


where e is an appropriate vector. Matrix C has two diagonal 
sub-matrices 

Jn » A H 


and 


C a = r'DAA' 


C b = DAY, 


and the other two sub-matrices can be represented as 


and 


C c = r'V A DAA 

C d = 51 - ri +1 \ H A ff DAV. 


(19) 

( 20 ) 

( 21 ) 

( 22 ) 


Letting uj = 0, we have 

N 

Kn) = lliv - ft)(ft - 5) + + -Vft 2 - (28) 


T t +1 r *+l 


2=1 


As the variance always converges, we have r* +1 = = t x 

after a certain number of iterations. Then d28t can be reduced 
to 


N 


Kn) = P[(?7 2 -aty + aft) 


i =1 


where from (fl7l) 


min{M,N} . 

_ l Tc|a 

~~ AT ^ 


TV ' T, 

2=1 


|Ai| 2 + cr 


2 ' 


(29) 


(30) 


which indicates that 0 < a < 1. Hence the eigenvalues are 
given by 


Vi( 1 , 2 ) — 


—a ± yja 2 — 4aft 


(31) 


for i = 1,2, ...,7V. Recall that 0 < ft < 1. If a > /3j/4, the 
eigenvalues are real, and it can be easily shown that 


Next, we find the eigenvalues of matrix C, i.e., the roots of 
the following polynomial 


h( V ) = \ V l-C\ 


rjl ~ C a C b 
C c rjl - C d 


(23) 


We note that the identity matrices in (|23| > have different sizes 
(i.e., the use of I is abused for notation simplification). As C a 
is a diagonal matrix (with non-negative elements), a diagonal 
matrix wl can be used to guarantee that //I — C„ + +1 is 
invertible. Define a new polynomial 


h a (v) 


r]I-C a + wl C b 

C c rjl - C d 


(24) 


Clearly the roots of h a (rj) with u> = 0 are the eigenvalues of 
matrix C. It can be shown that h a {rj) can be rewritten as 


h a { V )= |? ? I - C Q + wl| x |t,I -C d - C c {rjl -C a + wI)- 1 C b | 
=\rjl — r*DAA ff + wl| x 

\\ H \ x |(77 - 5)1 + r* +1 A ff DA + r* +1 r*A fl DAA ff 
( 77 I - r*DAA fl + wI) _1 DA| x |V|. (25) 


As V is a unitary matrix, |V H | = |V| = 1. So they can 
be removed from (l25l >. Note that A is a rectangular diagonal 
matrix with size M x TV and all the matrices left in (l25l ) are 
diagonal. 

Case 1: M = TV. In this case, A is a diagonal matrix. Define 
vector (3 = [/3i,..., (3n] T whose elements are the diagonal 
elements of r*DAA ff , i.e.. 


ft = 


ilftl 


(26) 


r*|Ai | 2 + CT 2 ' 

where we can see that 0 < fa < 1. Hence d25l) can be rewritten 
as 


N 

h a (v)=Y[(v - Pi 
2=1 
N 


=YI(V ~ Pi + w)((v 

2=1 


a) + 

5) + 


r t~\~ 1 


r t + l 


(ft + 

ft) + 


ft 2 


rj - Pi + w 


)) 


r i+l 


-Pi (27) 


\r]i\ < ol < 1, (32) 

where the equality holds when ft = 0. If a < Pi/A, the 
eigenvalues are complex valued, and it can be shown that 

\Vi\ < a Pi < 1- (33) 

Case 2: M > TV. In this case, A is a ‘tall’ rectangular 
diagonal matrix. We define diagonal matrix A with size TV x TV 
as the upper part of A, and define diagonal matrix D with size 
TV x TV as the upper left part of D (whose size is M x M). It 
is not hard to show that 

h a (r]) = \rjl M -N +ujI m ~n\ x ftl - t*DAA h + wl| 

x | (r) - a)I + r* +1 A h DA + r'+V'A^DAA® 
{rjl - t‘DAA h + wI) -1 DA|. (34) 

After some manipulations, we have 

N 

h{rj) = r/ M ~ N Win 2 -&n + aft). (35) 

2=1 

Hence the eigenvalues are the same as those in Case 1 except 
that M — TV eigenvalues are zero. 

Case 3: M < TV. In this case, A is a ‘fat’ rectangular 
diagonal matrix. Define diagonal matrix A with sizeM x M 
as the left part of A. We can show that 

h a {rj) = \rjl - t)DAA H + wl| x | (77 - a)I + r‘ +1 A ff DA + 
r* +1 r‘A"DAA H (77l - r*DAA H + wI)“ 1 DA| 
x | (77 — a)Ijv-M|- (36) 

Then, we can have 

M 

h{n) = {n~Oi) N ~ M ff(?7 2 -an + api). (37) 

2=1 

The eigenvalues are the same as those in Case 1 except that 
TV — M eigenvalues are a. 

The above shows that \n, \ < a for all the cases (noting that 
aft < a) and any matrix A. Because a is smaller than 1, the 
algorithm converges for any A. 
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